// Licensed to the .NET Foundation under one or more agreements.
// The .NET Foundation licenses this file to you under the MIT license.
// See the LICENSE file in the project root for more information.

using System;
using System.Diagnostics;
using System.Collections.Generic;
using Microsoft.ML.Probabilistic;
using Microsoft.ML.Probabilistic.Utilities;
using Microsoft.ML.Probabilistic.Math;
using Microsoft.ML.Probabilistic.Distributions;
using Microsoft.ML.Probabilistic.Collections;
using Microsoft.ML.Probabilistic.Factors;

/*
namespace Microsoft.ML.Probabilistic.Tests
{
    /// <summary>
    /// This is a test usage case of the nonconjugate messages used by the ProductExp factor. 
    /// </summary>
    /// <remarks>
    /// A simple test can be run using
    /// (new LogisticIRT_VMP()).Test(); 
    /// 
    /// The easiest way to use this class is to wrap an instance in a var object and use
    /// the methods on var to set parameters and execute inference.
    /// 
    /// If you instead wish to use this class directly, you must perform the following steps:
    /// 1) Create an instance of the class
    /// 2) Set the value of any externally-set fields e.g. data, priors
    /// 3) Call the Execute(numberOfIterations) method
    /// 4) Use the XXXMarginal() methods to retrieve posterior marginals for different variables.
    /// 
    /// Generated by Infer.NET 2.4 beta 1 at 12:08 on 22 March 2010.
    /// </remarks>
    public class LogisticIRT_VMP
    {
        #region Fields
        public int NumQuestions;
        public int NumStudents;
        /// <summary>Field backing the response property</summary>
        private bool[,] Response;
        /// <summary>The number of iterations last computed by Constant. Set this to zero to force re-execution of Constant</summary>
        public int Constant_iterationsDone;
        /// <summary>The number of iterations last computed by Init_response. Set this to zero to force re-execution of Init_response</summary>
        public int Init_response_iterationsDone;
        /// <summary>True if Init_response has performed initialisation. Set this to false to force re-execution of Init_response</summary>
        public bool Init_response_isInitialised;
        /// <summary>The number of iterations last computed by Changed_response. Set this to zero to force re-execution of Changed_response</summary>
        public int Changed_response_iterationsDone;
        /// <summary>Field backing the AfterUpdate property</summary>
        private Action<int> afterUpdate;
        /// <summary>Field backing the ResumeLastRun property</summary>
        private bool resumeLastRun;
        /// <summary>The constant 'vGaussian2'</summary>
        public Gaussian LogDiscPrior;
        /// <summary>Message from definition of 'logDisc'</summary>
        public DistributionStructArray<Gaussian,double> logDisc_F;
        /// <summary>Messages to uses of 'logDisc'</summary>
        public DistributionStructArray<Gaussian,double>[] logDisc_uses_F;
        /// <summary>Messages from uses of 'logDisc'</summary>
        public DistributionStructArray<NonconjugateGaussian,double>[] logDisc_uses_B;
        /// <summary>Messages from uses of 'minus'</summary>
        public DistributionStructArray2D<Gaussian,double>[] minus_uses_B;
        /// <summary>The constant 'vGaussian0'</summary>
        public Gaussian AbilityPrior;
        /// <summary>Message to definition of 'weighted'</summary>
        public DistributionStructArray2D<Gaussian,double> weighted_B;
        /// <summary>Message to definition of 'minus'</summary>
        public DistributionStructArray2D<Gaussian,double> minus_B;
        /// <summary>Message to definition of 'difficulty_rep0'</summary>
        public DistributionStructArray<Gaussian,double>[] difficulty_rep0_B;
        /// <summary>The constant 'vGaussian1'</summary>
        public Gaussian DifficultyPrior;
        /// <summary>Messages to uses of 'ability'</summary>
        public DistributionStructArray<Gaussian,double>[] ability_uses_F;
        /// <summary>Messages from uses of 'difficulty'</summary>
        public DistributionStructArray<Gaussian,double>[] difficulty_uses_B;
        /// <summary>Message from definition of 'difficulty'</summary>
        public DistributionStructArray<Gaussian,double> difficulty_F;
        /// <summary>Messages to uses of 'difficulty'</summary>
        public DistributionStructArray<Gaussian,double>[] difficulty_uses_F;
        /// <summary>Message to definition of 'ability_rep0'</summary>
        public DistributionStructArray<Gaussian,double>[] ability_rep0_B;
        /// <summary>Messages from uses of 'ability'</summary>
        public DistributionStructArray<Gaussian,double>[] ability_uses_B;
        /// <summary>Message from definition of 'ability'</summary>
        public DistributionStructArray<Gaussian,double> ability_F;
        /// <summary>Message from definition of 'minus'</summary>
        public DistributionStructArray2D<Gaussian,double> minus_F;
        /// <summary>Message from definition of 'weighted'</summary>
        public DistributionStructArray2D<Gaussian,double> weighted_F;
        /// <summary>Messages from uses of 'weighted'</summary>
        public DistributionStructArray2D<Gaussian,double>[] weighted_uses_B;
        /// <summary>Message to definition of 'logDisc_rep0'</summary>
        public DistributionStructArray<NonconjugateGaussian,double>[] logDisc_rep0_B;
        /// <summary>Message to marginal of 'logDisc'</summary>
        public DistributionStructArray<Gaussian,double> logDisc_marginal_B;
        /// <summary>Message to marginal of 'ability'</summary>
        public DistributionStructArray<Gaussian,double> ability_marginal_B;
        /// <summary>Message to marginal of 'difficulty'</summary>
        public DistributionStructArray<Gaussian,double> difficulty_marginal_B;
        /// <summary>Message to marginal of 'minus'</summary>
        public DistributionStructArray2D<Gaussian,double> minus_marginal_B;
        /// <summary>Message to marginal of 'weighted'</summary>
        public DistributionStructArray2D<Gaussian,double> weighted_marginal_B;
        #endregion

        public void Test()
        {
            //NumStudents = 10;
            //NumQuestions = 12;
            //bool[,] data = new bool[NumStudents, NumQuestions];
            //for (int i = 0; i < NumStudents; i++)
            //{
            //    for (int j = 0; j < NumQuestions; j++)
            //    {
            //        data[i, j] = (i > j);
            //    }
            //}
            //response = data;

            Dictionary<string, object> dict = MatlabReader.Read(@"..\..\..\Tests\Data\IRT2PL_10_1000.mat");
            Matrix m = (Matrix)dict["Y"];
            var data = PsychTests.ConvertToBool(m.ToArray());
            m = (Matrix)dict["discrimination"];
            var discriminationTrue = Util.ArrayInit(data.GetLength(1), i => m[i]);
            m = (Matrix)dict["ability"];
            var abilityTrue = Util.ArrayInit(data.GetLength(0), i => m[i]);
            m = (Matrix)dict["difficulty"];
            var difficultyTrue = Util.ArrayInit(data.GetLength(1), i => m[i]);
            NumStudents = data.GetLength(0);
            NumQuestions = data.GetLength(1);
            response = data;
            var s = new Stopwatch(); 
            s.Start();
            Execute(50);
            s.Stop(); 

            var abilityMarg = AbilityMarginal(); 
            Console.WriteLine("Abilities:");
            Console.WriteLine("Inferred\tTruth");
            for (int i = 0; i < NumStudents; i++)
                Console.WriteLine(abilityMarg[i].GetMean() + "\t" + abilityTrue[i]);

            var difficultyMarg = DifficultyMarginal(); 
            Console.WriteLine("Difficulties:");
            Console.WriteLine("Inferred\tTruth");
            for (int i = 0; i < NumQuestions; i++)
                Console.WriteLine(difficultyMarg[i].GetMean() + "\t" + difficultyTrue[i]);

            Console.WriteLine("Compare inferred logDiscrimination to ground truth:"); 
            var marg = LogDiscMarginal();
            for (int i=0;i<NumQuestions;i++)
                Console.WriteLine(marg[i].GetMean() + " should be " + Math.Log(discriminationTrue[i]));

            Console.WriteLine("Inference took {0} ms", s.ElapsedMilliseconds); 
        }

        #region Properties
        /// <summary>The externally-specified value of 'response'</summary>
        public bool[,] response
{            get {
                return this.Response;
            }
            set {
                if (this.Response!=value) {
                    this.Response = value;
                    this.Init_response_isInitialised = false;
                    this.Changed_response_iterationsDone = 0;
                }
            }
        }

        /// <summary>Called after each iteration</summary>
        public Action<int> AfterUpdate
{            get {
                return this.afterUpdate;
            }
            set {
                this.afterUpdate = value;
            }
        }

        /// <summary>Set to true to re-use previous message state</summary>
        public bool ResumeLastRun
{            get {
                return this.resumeLastRun;
            }
            set {
                this.resumeLastRun = value;
            }
        }

        #endregion

        #region Methods
        /// <summary>Update all marginals, by iterating message passing the given number of times</summary>
        /// <param name="numberOfIterations">The number of times to iterate each loop</param>
        public void Execute(int numberOfIterations)
        {
            this.Constant();
            this.Init_response();
            this.Changed_response(numberOfIterations);
        }

        /// <summary>Reset all messages to their initial values on the next call to Execute</summary>
        public void Reset()
        {
            this.Constant_iterationsDone = 0;
        }

        /// <summary>Computations that do not depend on observed values</summary>
        public void Constant()
        {
            if (this.Constant_iterationsDone==1) {
                return ;
            }
            this.LogDiscPrior = Gaussian.FromNatural(0, 1.5e-2);
            //this.NumQuestions = 2;
            // Create array for 'logDisc' forwards messages.
            this.logDisc_F = new DistributionStructArray<Gaussian,double>(this.NumQuestions);
            for(int questions = 0; questions<NumQuestions; questions++) {
                this.logDisc_F[questions] = ArrayHelper.MakeUniform<Gaussian>(this.LogDiscPrior);
                // Message to 'logDisc' from Random factor
                this.logDisc_F[questions] = UnaryOp<double>.RandomAverageLogarithm<Gaussian>(this.LogDiscPrior);
            }
            // Create array for 'logDisc_uses' forwards messages.
            this.logDisc_uses_F = new DistributionStructArray<Gaussian,double>[1];
            for(int _ind = 0; _ind<1; _ind++) {
                // Create array for 'logDisc_uses' forwards messages.
                this.logDisc_uses_F[_ind] = new DistributionStructArray<Gaussian,double>(this.NumQuestions);
                for(int questions = 0; questions<NumQuestions; questions++) {
                    this.logDisc_uses_F[_ind][questions] = ArrayHelper.MakeUniform<Gaussian>(this.LogDiscPrior);
                }
            }
            // Create array for 'logDisc_uses' backwards messages.
            this.logDisc_uses_B = new DistributionStructArray<NonconjugateGaussian,double>[1];
            for(int _ind = 0; _ind<1; _ind++) {
                // Create array for 'logDisc_uses' backwards messages.
                this.logDisc_uses_B[_ind] = new DistributionStructArray<NonconjugateGaussian, double>(this.NumQuestions);
            }
            //this.NumStudents = 20;
            // Create array for 'minus_uses' backwards messages.
            this.minus_uses_B = new DistributionStructArray2D<Gaussian,double>[1];
            for(int _ind = 0; _ind<1; _ind++) {
                // Create array for 'minus_uses' backwards messages.
                this.minus_uses_B[_ind] = new DistributionStructArray2D<Gaussian,double>(this.NumStudents, this.NumQuestions);
            }
            this.AbilityPrior = Gaussian.FromNatural(0, 1E-6);
            for(int _ind = 0; _ind<1; _ind++) {
                for(int students = 0; students<NumStudents; students++) {
                    for(int questions = 0; questions<NumQuestions; questions++) {
                        this.minus_uses_B[_ind][students, questions] = ArrayHelper.MakeUniform<Gaussian>(this.AbilityPrior);
                    }
                }
            }
            // Create array for 'weighted' backwards messages.
            this.weighted_B = new DistributionStructArray2D<Gaussian,double>(this.NumStudents, this.NumQuestions);
            for(int students = 0; students<NumStudents; students++) {
                for(int questions = 0; questions<NumQuestions; questions++) {
                    this.weighted_B[students, questions] = ArrayHelper.MakeUniform<Gaussian>(new Gaussian());
                }
            }
            // Create array for 'minus' backwards messages.
            this.minus_B = new DistributionStructArray2D<Gaussian,double>(this.NumStudents, this.NumQuestions);
            for(int students = 0; students<NumStudents; students++) {
                for(int questions = 0; questions<NumQuestions; questions++) {
                    this.minus_B[students, questions] = ArrayHelper.MakeUniform<Gaussian>(this.AbilityPrior);
                }
            }
            // Create array for replicates of 'difficulty_rep0_B'
            this.difficulty_rep0_B = new DistributionStructArray<Gaussian,double>[NumQuestions];
            for(int questions = 0; questions<NumQuestions; questions++) {
                // Create array for 'difficulty_rep0' backwards messages.
                this.difficulty_rep0_B[questions] = new DistributionStructArray<Gaussian,double>(this.NumStudents);
            }
            this.DifficultyPrior = Gaussian.FromNatural(0, 1E-6);
            for(int questions = 0; questions<NumQuestions; questions++) {
                for(int students = 0; students<NumStudents; students++) {
                    this.difficulty_rep0_B[questions][students] = ArrayHelper.MakeUniform<Gaussian>(this.DifficultyPrior);
                }
            }
            // Create array for 'ability_uses' forwards messages.
            this.ability_uses_F = new DistributionStructArray<Gaussian,double>[1];
            for(int _ind = 0; _ind<1; _ind++) {
                // Create array for 'ability_uses' forwards messages.
                this.ability_uses_F[_ind] = new DistributionStructArray<Gaussian,double>(this.NumStudents);
                for(int students = 0; students<NumStudents; students++) {
                    this.ability_uses_F[_ind][students] = ArrayHelper.MakeUniform<Gaussian>(this.AbilityPrior);
                }
            }
            // Create array for 'difficulty_uses' backwards messages.
            this.difficulty_uses_B = new DistributionStructArray<Gaussian,double>[1];
            for(int _ind = 0; _ind<1; _ind++) {
                // Create array for 'difficulty_uses' backwards messages.
                this.difficulty_uses_B[_ind] = new DistributionStructArray<Gaussian,double>(this.NumQuestions);
            }
            // Create array for 'difficulty' forwards messages.
            this.difficulty_F = new DistributionStructArray<Gaussian,double>(this.NumQuestions);
            for(int questions = 0; questions<NumQuestions; questions++) {
                this.difficulty_F[questions] = ArrayHelper.MakeUniform<Gaussian>(this.DifficultyPrior);
                // Message to 'difficulty' from Random factor
                this.difficulty_F[questions] = UnaryOp<double>.RandomAverageLogarithm<Gaussian>(this.DifficultyPrior);
            }
            // Create array for 'difficulty_uses' forwards messages.
            this.difficulty_uses_F = new DistributionStructArray<Gaussian,double>[1];
            for(int _ind = 0; _ind<1; _ind++) {
                // Create array for 'difficulty_uses' forwards messages.
                this.difficulty_uses_F[_ind] = new DistributionStructArray<Gaussian,double>(this.NumQuestions);
                for(int questions = 0; questions<NumQuestions; questions++) {
                    this.difficulty_uses_F[_ind][questions] = ArrayHelper.MakeUniform<Gaussian>(this.DifficultyPrior);
                }
            }
            // Create array for replicates of 'ability_rep0_B'
            this.ability_rep0_B = new DistributionStructArray<Gaussian,double>[NumStudents];
            for(int students = 0; students<NumStudents; students++) {
                // Create array for 'ability_rep0' backwards messages.
                this.ability_rep0_B[students] = new DistributionStructArray<Gaussian,double>(this.NumQuestions);
                for(int questions = 0; questions<NumQuestions; questions++) {
                    this.ability_rep0_B[students][questions] = ArrayHelper.MakeUniform<Gaussian>(this.AbilityPrior);
                }
            }
            // Create array for 'ability_uses' backwards messages.
            this.ability_uses_B = new DistributionStructArray<Gaussian,double>[1];
            for(int _ind = 0; _ind<1; _ind++) {
                // Create array for 'ability_uses' backwards messages.
                this.ability_uses_B[_ind] = new DistributionStructArray<Gaussian,double>(this.NumStudents);
            }
            // Create array for 'ability' forwards messages.
            this.ability_F = new DistributionStructArray<Gaussian,double>(this.NumStudents);
            for(int students = 0; students<NumStudents; students++) {
                this.ability_F[students] = ArrayHelper.MakeUniform<Gaussian>(this.AbilityPrior);
                // Message to 'ability' from Random factor
                this.ability_F[students] = UnaryOp<double>.RandomAverageLogarithm<Gaussian>(this.AbilityPrior);
            }
            // Create array for 'minus' forwards messages.
            this.minus_F = new DistributionStructArray2D<Gaussian,double>(this.NumStudents, this.NumQuestions);
            for(int students = 0; students<NumStudents; students++) {
                for(int questions = 0; questions<NumQuestions; questions++) {
                    this.minus_F[students, questions] = ArrayHelper.MakeUniform<Gaussian>(this.AbilityPrior);
                }
            }
            // Create array for 'weighted' forwards messages.
            this.weighted_F = new DistributionStructArray2D<Gaussian,double>(this.NumStudents, this.NumQuestions);
            for(int students = 0; students<NumStudents; students++) {
                for(int questions = 0; questions<NumQuestions; questions++) {
                    this.weighted_F[students, questions] = ArrayHelper.MakeUniform<Gaussian>(new Gaussian());
                }
            }
            // Create array for 'weighted_uses' backwards messages.
            this.weighted_uses_B = new DistributionStructArray2D<Gaussian,double>[1];
            for(int _ind = 0; _ind<1; _ind++) {
                // Create array for 'weighted_uses' backwards messages.
                this.weighted_uses_B[_ind] = new DistributionStructArray2D<Gaussian,double>(this.NumStudents, this.NumQuestions);
                for(int students = 0; students<NumStudents; students++) {
                    for(int questions = 0; questions<NumQuestions; questions++) {
                        this.weighted_uses_B[_ind][students, questions] = ArrayHelper.MakeUniform<Gaussian>(new Gaussian());
                    }
                }
            }
            // Create array for replicates of 'logDisc_rep0_B'
            this.logDisc_rep0_B = new DistributionStructArray<NonconjugateGaussian, double>[NumQuestions];
            for(int questions = 0; questions<NumQuestions; questions++) {
                // Create array for 'logDisc_rep0' backwards messages.
                this.logDisc_rep0_B[questions] = new DistributionStructArray<NonconjugateGaussian, double>(this.NumStudents);
                for(int students = 0; students<NumStudents; students++) {
                    this.logDisc_rep0_B[questions][students] = ArrayHelper.MakeUniform<NonconjugateGaussian>(new NonconjugateGaussian());
                }
            }
            this.Constant_iterationsDone = 1;
            this.Init_response_iterationsDone = 0;
            this.Changed_response_iterationsDone = 0;
        }

        /// <summary>Computations that depend on the observed value of </summary>
        public void Init_response()
        {
            if ((this.Init_response_iterationsDone==1)&&(this.ResumeLastRun||this.Init_response_isInitialised)) {
                return ;
            }
            for(int _ind = 0; _ind<1; _ind++) {
                for(int questions = 0; questions<NumQuestions; questions++) {
                    this.logDisc_uses_B[_ind][questions] = ArrayHelper.MakeUniform<NonconjugateGaussian>(new NonconjugateGaussian());
                }
            }
            for(int _ind = 0; _ind<1; _ind++) {
                for(int questions = 0; questions<NumQuestions; questions++) {
                    this.difficulty_uses_B[_ind][questions] = ArrayHelper.MakeUniform<Gaussian>(this.DifficultyPrior);
                }
            }
            for(int _ind = 0; _ind<1; _ind++) {
                for(int students = 0; students<NumStudents; students++) {
                    this.ability_uses_B[_ind][students] = ArrayHelper.MakeUniform<Gaussian>(this.AbilityPrior);
                }
            }
            // Message to 'logDisc_uses' from UsesEqualDef factor
            this.logDisc_uses_F[0] = NonconjugateUsesEqualDefOp.UsesAverageLogarithm(this.logDisc_uses_B, this.logDisc_F, this.logDisc_uses_F[0]);
            // Message to 'difficulty_uses' from UsesEqualDef factor
            this.difficulty_uses_F[0] = UsesEqualDefVmpOp.UsesAverageLogarithm<DistributionStructArray<Gaussian,double>>(this.difficulty_uses_B, this.difficulty_F, 0, this.difficulty_uses_F[0]);
            // Message to 'ability_uses' from UsesEqualDef factor
            this.ability_uses_F[0] = UsesEqualDefVmpOp.UsesAverageLogarithm<DistributionStructArray<Gaussian,double>>(this.ability_uses_B, this.ability_F, 0, this.ability_uses_F[0]);
            for(int questions = 0; questions<NumQuestions; questions++) {
                for(int students = 0; students<NumStudents; students++) {
                    // Message to 'minus' from Difference factor
                    this.minus_F[students, questions] = DoubleMinusVmpOp.DifferenceAverageLogarithm(this.ability_uses_F[0][students], this.difficulty_uses_F[0][questions]);
                    // Message to 'weighted' from ProductExp factor
                    this.weighted_F[students, questions] = ProductExpOp.ProductExpAverageLogarithm(this.minus_F[students, questions], this.logDisc_uses_F[0][questions]);
                }
            }
            this.Init_response_iterationsDone = 1;
            this.Init_response_isInitialised = true;
            this.Changed_response_iterationsDone = 0;
        }

        /// <summary>Computations that depend on the observed value of response</summary>
        /// <param name="numberOfIterations">The number of times to iterate each loop</param>
        public void Changed_response(int numberOfIterations)
        {
            if (this.Changed_response_iterationsDone==numberOfIterations) {
                return ;
            }
            for(int iteration = 0; iteration<numberOfIterations; iteration++) {
                for(int questions = 0; questions<NumQuestions; questions++) {
                    for(int students = 0; students<NumStudents; students++) {
                        // Message to 'weighted_uses' from BernoulliFromLogOdds factor
                        this.weighted_uses_B[0][students, questions] = BernoulliFromLogOddsOp.LogOddsAverageLogarithm(this.response[students, questions], this.weighted_F[students, questions], this.weighted_uses_B[0][students, questions]);
                    }
                }
                // Message to 'weighted' from ReplicateWithMarginal factor
                this.weighted_B = ReplicateOp.DefAverageLogarithm<DistributionStructArray2D<Gaussian,double>>(this.weighted_uses_B, this.weighted_B);
                for(int questions = 0; questions<NumQuestions; questions++) {
                    for(int students = 0; students<NumStudents; students++) {
                        // Message to 'minus_uses' from ProductExp factor
                        this.minus_uses_B[0][students, questions] = ProductExpOp.AAverageLogarithm(this.weighted_B[students, questions], this.logDisc_uses_F[0][questions]);
                    }
                }
                // Message to 'minus' from ReplicateWithMarginal factor
                this.minus_B = ReplicateOp.DefAverageLogarithm<DistributionStructArray2D<Gaussian,double>>(this.minus_uses_B, this.minus_B);
                for(int questions = 0; questions<NumQuestions; questions++) {
                    for(int students = 0; students<NumStudents; students++) {
                        // Message to 'ability_rep0' from Difference factor
                        this.ability_rep0_B[students][questions] = DoubleMinusVmpOp.AAverageLogarithm(this.minus_B[students, questions], this.difficulty_uses_F[0][questions]);
                    }
                }
                for(int students = 0; students<NumStudents; students++) {
                    // Message to 'ability_uses' from Replicate factor
                    this.ability_uses_B[0][students] = ReplicateOp.DefAverageLogarithm<Gaussian>(this.ability_rep0_B[students], this.ability_uses_B[0][students]);
                }
                // Message to 'ability_uses' from UsesEqualDef factor
                this.ability_uses_F[0] = UsesEqualDefVmpOp.UsesAverageLogarithm<DistributionStructArray<Gaussian,double>>(this.ability_uses_B, this.ability_F, 0, this.ability_uses_F[0]);
                for(int questions = 0; questions<NumQuestions; questions++) {
                    for(int students = 0; students<NumStudents; students++) {
                        // Message to 'difficulty_rep0' from Difference factor
                        this.difficulty_rep0_B[questions][students] = DoubleMinusVmpOp.BAverageLogarithm(this.minus_B[students, questions], this.ability_uses_F[0][students]);
                    }
                    // Message to 'difficulty_uses' from Replicate factor
                    this.difficulty_uses_B[0][questions] = ReplicateOp.DefAverageLogarithm<Gaussian>(this.difficulty_rep0_B[questions], this.difficulty_uses_B[0][questions]);
                }
                // Message to 'difficulty_uses' from UsesEqualDef factor
                this.difficulty_uses_F[0] = UsesEqualDefVmpOp.UsesAverageLogarithm<DistributionStructArray<Gaussian,double>>(this.difficulty_uses_B, this.difficulty_F, 0, this.difficulty_uses_F[0]);
                for(int questions = 0; questions<NumQuestions; questions++) {
                    for(int students = 0; students<NumStudents; students++) {
                        // Message to 'minus' from Difference factor
                        this.minus_F[students, questions] = DoubleMinusVmpOp.DifferenceAverageLogarithm(this.ability_uses_F[0][students], this.difficulty_uses_F[0][questions]);
                        // Message to 'logDisc_rep0' from ProductExp factor
                        this.logDisc_rep0_B[questions][students] = ProductExpOp.BAverageLogarithm(this.weighted_B[students, questions], this.minus_F[students, questions], this.logDisc_uses_F[0][questions], this.logDisc_rep0_B[questions][students]);
                    }
                    // Message to 'logDisc_uses' from Replicate factor
                    this.logDisc_uses_B[0][questions] = ReplicateOp.DefAverageLogarithm<NonconjugateGaussian>(this.logDisc_rep0_B[questions], this.logDisc_uses_B[0][questions]);
                }
                // Message to 'logDisc_uses' from UsesEqualDef factor
                this.logDisc_uses_F[0] = NonconjugateUsesEqualDefOp.UsesAverageLogarithm(this.logDisc_uses_B, this.logDisc_F, this.logDisc_uses_F[0]);
                for(int questions = 0; questions<NumQuestions; questions++) {
                    for(int students = 0; students<NumStudents; students++) {
                        // Message to 'weighted' from ProductExp factor
                        this.weighted_F[students, questions] = ProductExpOp.ProductExpAverageLogarithm(this.minus_F[students, questions], this.logDisc_uses_F[0][questions]);
                    }
                }
                if (this.AfterUpdate!=null) {
                    this.AfterUpdate.Invoke(iteration);
                }
            }
            // Create array for 'logDisc_marginal' backwards messages.
            this.logDisc_marginal_B = new DistributionStructArray<Gaussian,double>(this.NumQuestions);
            for(int questions = 0; questions<NumQuestions; questions++) {
                this.logDisc_marginal_B[questions] = ArrayHelper.MakeUniform<Gaussian>(this.LogDiscPrior);
            }
            // Message to 'logDisc_marginal' from UsesEqualDef factor
            this.logDisc_marginal_B = NonconjugateUsesEqualDefOp.MarginalAverageLogarithm(this.logDisc_uses_B, this.logDisc_F, this.logDisc_marginal_B);
            // Create array for 'ability_marginal' backwards messages.
            this.ability_marginal_B = new DistributionStructArray<Gaussian,double>(this.NumStudents);
            for(int students = 0; students<NumStudents; students++) {
                this.ability_marginal_B[students] = ArrayHelper.MakeUniform<Gaussian>(this.AbilityPrior);
            }
            // Message to 'ability_marginal' from UsesEqualDef factor
            this.ability_marginal_B = UsesEqualDefVmpOp.MarginalAverageLogarithm<DistributionStructArray<Gaussian,double>>(this.ability_uses_B, this.ability_F, this.ability_marginal_B);
            // Create array for 'difficulty_marginal' backwards messages.
            this.difficulty_marginal_B = new DistributionStructArray<Gaussian,double>(this.NumQuestions);
            for(int questions = 0; questions<NumQuestions; questions++) {
                this.difficulty_marginal_B[questions] = ArrayHelper.MakeUniform<Gaussian>(this.DifficultyPrior);
            }
            // Message to 'difficulty_marginal' from UsesEqualDef factor
            this.difficulty_marginal_B = UsesEqualDefVmpOp.MarginalAverageLogarithm<DistributionStructArray<Gaussian,double>>(this.difficulty_uses_B, this.difficulty_F, this.difficulty_marginal_B);
            // Create array for 'minus_marginal' backwards messages.
            this.minus_marginal_B = new DistributionStructArray2D<Gaussian,double>(this.NumStudents, this.NumQuestions);
            for(int students = 0; students<NumStudents; students++) {
                for(int questions = 0; questions<NumQuestions; questions++) {
                    this.minus_marginal_B[students, questions] = ArrayHelper.MakeUniform<Gaussian>(this.AbilityPrior);
                }
            }
            // Message to 'minus_marginal' from ReplicateWithMarginal factor
            this.minus_marginal_B = ReplicateOp.MarginalAverageLogarithm<DistributionStructArray2D<Gaussian,double>,DistributionStructArray2D<Gaussian,double>>(this.minus_F, this.minus_marginal_B);
            // Create array for 'weighted_marginal' backwards messages.
            this.weighted_marginal_B = new DistributionStructArray2D<Gaussian,double>(this.NumStudents, this.NumQuestions);
            for(int students = 0; students<NumStudents; students++) {
                for(int questions = 0; questions<NumQuestions; questions++) {
                    this.weighted_marginal_B[students, questions] = ArrayHelper.MakeUniform<Gaussian>(new Gaussian());
                }
            }
            // Message to 'weighted_marginal' from ReplicateWithMarginal factor
            this.weighted_marginal_B = ReplicateOp.MarginalAverageLogarithm<DistributionStructArray2D<Gaussian,double>,DistributionStructArray2D<Gaussian,double>>(this.weighted_F, this.weighted_marginal_B);
            this.Changed_response_iterationsDone = numberOfIterations;
        }

        /// <summary>
        /// Returns the marginal distribution for 'logDisc' given by the current state of the
        /// message passing algorithm.
        /// </summary>
        /// <returns>The marginal distribution</returns>
        public DistributionStructArray<Gaussian,double> LogDiscMarginal()
        {
            return this.logDisc_marginal_B;
        }

        /// <summary>
        /// Returns the marginal distribution for 'ability' given by the current state of the
        /// message passing algorithm.
        /// </summary>
        /// <returns>The marginal distribution</returns>
        public DistributionStructArray<Gaussian,double> AbilityMarginal()
        {
            return this.ability_marginal_B;
        }

        /// <summary>
        /// Returns the marginal distribution for 'difficulty' given by the current state of the
        /// message passing algorithm.
        /// </summary>
        /// <returns>The marginal distribution</returns>
        public DistributionStructArray<Gaussian,double> DifficultyMarginal()
        {
            return this.difficulty_marginal_B;
        }

        /// <summary>
        /// Returns the marginal distribution for 'minus' given by the current state of the
        /// message passing algorithm.
        /// </summary>
        /// <returns>The marginal distribution</returns>
        public DistributionStructArray2D<Gaussian,double> MinusMarginal()
        {
            return this.minus_marginal_B;
        }

        /// <summary>
        /// Returns the marginal distribution for 'weighted' given by the current state of the
        /// message passing algorithm.
        /// </summary>
        /// <returns>The marginal distribution</returns>
        public DistributionStructArray2D<Gaussian,double> WeightedMarginal()
        {
            return this.weighted_marginal_B;
        }

        #endregion

    }
}
*/